subroutine CalcCorrectorStep()
use Data_Cfoam, only: RhoNew,Unew,Pnew,Enew, Rho,U,E,P, cfl
use Data_Mesh, only: nCells, meshV, DbCells, meshCells, meshFaceOwner, nInternalFaces
use DataCells, only: sumCellFluxes, fluxesCellPresent, SubStepCellVars, CharsCellPresent, &
                     CellGlobalStepDone, CurrentVelocities
use DataTime, only: TimeCell, dtCell, dt2Cell, dtFace, dtGlobal=>dtOut, &
                    TimeGlobal=>ttime, &
                    TimeStepMin,TimeStepMax, TimeStepDiff, &
                    NextTimeCell, CurrentTimeCell, Balans0, TimeEpsilon, &
                    nNextTimeCell,nNextTimeLocal, nTimeStepDiff, &
                    AdaptiveMeshRefinement, nDtGlobal, dnTimeCell
use DataPerfectGas, only: PressureAllVars, SoundVelocityPRho
use DataFaces, only: CharsFacePresent, fluxesFacePresent, CharsNew,CharsOld, fluxFaces, &
                     FaceVarsOld,FaceVarsNew, fluxFacesReserve, fluxCorrection, ViscousFluxes, &
                     FluxConserv
use DataViscosity, only: CalcViscosityOn
use DataSource
implicit none
    
    integer :: iCell, numFace, iFace, iCellOwner, iSide, n, i, m
    real(8) :: FluxesDivMeshV(5), Umod,Sound
    real(8) :: dtmin
    real(8) :: DivVStress, dtdV

    real(8) :: Pres, SumFluxCorrection(5), sumViscousCellFluxes(3)
    
    real(8) :: Pi,ek,a,Q, Rho0

    do iCell = 1,nCells
    
        if(CellGlobalStepDone(iCell)) cycle

        fluxesCellPresent(iCell) = .true.
        do numFace = 1,6
            iFace = meshCells(numFace,iCell)
            iCellOwner = meshFaceOwner(iFace)
            if(iCellOwner.eq.iCell)then
                iSide = 1
            else
                iSide = 2
            end if
            if(fluxesFacePresent(iSide,iFace)) cycle

            fluxesCellPresent(iCell) = .false.
            exit
        end do

        if(fluxesCellPresent(iCell).eqv..false.) cycle

        TimeCell(iCell) = TimeCell(iCell) + dtCell(iCell)

!        call CalcViscousFluxesForCell(iCell)

        sumCellFluxes(1:5,iCell) = 0d0
        SumFluxCorrection(1:5) = 0d0
        sumViscousCellFluxes(1:3) = 0d0
        DivVStress = 0.0d0
        do numFace = 1,6
            iFace = meshCells(numFace,iCell)
            iCellOwner = meshFaceOwner(iFace)
            if(iCellOwner.eq.iCell)then
                iSide = 1
            else
                iSide = 2
            end if

!if(iFace.le.nInternalFaces)then
!    FluxConserv(1:5,iSide,iFace) = FluxConserv(1:5,iSide,iFace) + &
!            dt2Cell(iCell) * fluxFaces(1:5,iSide,iFace) + fluxCorrection(1:5,iSide,iFace)
!end if

            sumCellFluxes(1:5,iCell) = sumCellFluxes(1:5,iCell) + fluxFaces(1:5,iSide,iFace)
            sumViscousCellFluxes(1:3) = sumViscousCellFluxes(1:3) + ViscousFluxes(1:3,iSide,iFace)
            SumFluxCorrection(1:5) = SumFluxCorrection(1:5) + fluxCorrection(1:5,iSide,iFace)

			do m = 1,3
    			DivVStress = DivVStress + FaceVarsNew(m+1,iSide,iFace) * ViscousFluxes(m,iSide,iFace)
			end do

            fluxCorrection(1:5,iSide,iFace) = 0d0

!            CharsFacePresent(iSide,iFace) = .false.
!            if(TimeGlobal-TimeCell(iCell).gt.TimeEpsilon)then
            if(dtGlobal-TimeCell(iCell).gt.TimeEpsilon)then
                fluxesFacePresent(iSide,iFace) = .false.
            end if
!            if(iFace.gt.nInternalFaces)then
!                Balans0(1:5) = Balans0(1:5) + fluxFaces(1:5,iSide,iFace) * dt2Cell(iCell)
!            end if
        end do

!        if(TimeGlobal-TimeCell(iCell).gt.TimeEpsilon)then
        if(dtGlobal-TimeCell(iCell).gt.TimeEpsilon)then
            CharsCellPresent(iCell) = .false.
        end if

        dtdV = dt2Cell(iCell) / meshV(iCell)

		SubStepCellVars(2:5,iCell) = SubStepCellVars(1,iCell) * SubStepCellVars(2:5,iCell) - &
                                     (sumCellFluxes(2:5,iCell) * dt2Cell(iCell) + SumFluxCorrection(2:5)) &
                                     / meshV(iCell)
		SubStepCellVars(2:4,iCell) = SubStepCellVars(2:4,iCell) + sumViscousCellFluxes(1:3) * dtdV ! dt2Cell(iCell) / meshV(iCell)
        SubStepCellVars(1,iCell) = SubStepCellVars(1,iCell) - (sumCellFluxes(1,iCell) * dt2Cell(iCell) + &
                                   SumFluxCorrection(1)) / meshV(iCell)
        SubStepCellVars(5,iCell) = SubStepCellVars(5,iCell) + DivVStress * dtdV

!	    if(SourcePresent)then
!	        if(iCell.eq.iCellSource)then
!		        call CalcSource(TimeCell(iCell)+TimeGlobal,Q) ! + dt2Cell(iCell)
!		        Q = Q * dt2Cell(iCell) / meshV(iCell)
!                SubStepCellVars(1,iCell) = SubStepCellVars(1,iCell) + Q
!            end if
!        end if

        SubStepCellVars(2:5,iCell) = SubStepCellVars(2:5,iCell) / SubStepCellVars(1,iCell)

!	    if(SourcePresent)then
!	        if(iCell.eq.iCellSource)then
!	            Rho0 = SubStepCellVars(1,iCell)
!!	            Q = -Q
!	            do i = 2,5
!	                SubStepCellVars(i,iCell) = SubStepCellVars(i,iCell) * Rho0 / (Rho0 - Q)
!                end do
!!	            SubStepCellVars(2,iCell) = SubStepCellVars(2,iCell) * Rho0 + Q * Unew(1,iCell)
!!		        SubStepCellVars(3,iCell) = SubStepCellVars(3,iCell) * Rho0 + Q * Unew(2,iCell)
!!		        SubStepCellVars(4,iCell) = SubStepCellVars(4,iCell) * Rho0 + Q * Unew(3,iCell)
!!		        SubStepCellVars(5,iCell) = SubStepCellVars(5,iCell) * Rho0 + Q * Enew(iCell)
!
!!                write(18,'(11f20.10)') TimeCell(iCell)+TimeGlobal,Q, Q*SubStepCellVars(2:5,iCell),SubStepCellVars(1:5,iCell) !  + dt2Cell(iCell)
!            end if
!        end if
!        
        if(SourcePresent)then
	        if(iCell.eq.iCellSource)then
		        call CalcSource(TimeCell(iCell)+TimeGlobal,Q) ! + dt2Cell(iCell)
		        Q = Q * dt2Cell(iCell) / meshV(iCell)
                SubStepCellVars(5,iCell) = SubStepCellVars(5,iCell) + Q / SubStepCellVars(1,iCell)
            end if
        end if
	        

!        CurrentTimeCell(iCell) = CurrentTimeCell(iCell) + dt2Cell(iCell)

!        TimeCell(iCell) = TimeCell(iCell) + dt2Cell(iCell)
        
!!		if(iCell.eq.1225)then ! 1760)then ! 1875)then
!		if(iCell.eq.28296)then ! 6957)then ! .or.iCell.eq.6958.or.iCell.eq.7117.or.iCell.eq.7118)then ! 1760)then ! 1875)then
!			Pi = 3.1415926535897932384626433832795
!		    ek = 1d7 * 1.739382343d0 / meshV(iCell)
!			a = 10.0d0 
!			Q = ek * sin(2.0d0 * a * Pi * TimeCell(iCell))
!			SubStepCellVars(5,iCell) = SubStepCellVars(5,iCell) + Q*dt2Cell(iCell)/SubStepCellVars(1,iCell)
!!			Enew(i) = Enew(i) + Q*dt2/rhonew(i)
!		end if

        ! 1:5 = Rho, Velocities, Energy
        Rho(iCell) = RhoNew(iCell)
        U(1:3,iCell) = Unew(1:3,iCell)
        E(iCell) = Enew(iCell)
        P(iCell) = Pnew(iCell)

        RhoNew(iCell) = SubStepCellVars(1,iCell)
        Unew(1:3,iCell) = SubStepCellVars(2:4,iCell)
        Enew(iCell) = SubStepCellVars(5,iCell)
        Pnew(iCell) = PressureAllVars(SubStepCellVars(1:5,iCell))
        
        U(1:3,iCell) = Unew(1:3,iCell)

        if(CalcViscosityOn)then
            do i =1,3
                CurrentVelocities(i,iCell) = Unew(i,iCell)
            end do
	    end if

        Umod = dsqrt(sum(Unew(1:3,iCell)*Unew(1:3,iCell)))
        Sound = SoundVelocityPRho(Pnew(iCell),RhoNew(iCell))

		dtCell(iCell) = DbCells(iCell) / (Umod + Sound)

        do numFace = 1,6
            iFace = meshCells(numFace,iCell)
            dtCell(iCell) = min(dtCell(iCell),dtFace(iFace))
            iCellOwner = meshFaceOwner(iFace)
            if(iCellOwner.eq.iCell)then
                iSide = 1
            else
                iSide = 2
            end if
            CharsOld(1:5,iSide,iFace) = CharsNew(1:5,iSide,iFace)
            FaceVarsOld(1:5,iSide,iFace) = FaceVarsNew(1:5,iSide,iFace) ! New2Old
        end do

        dtCell(iCell) = cfl * dtCell(iCell)
        dtCell(iCell) = min(dtCell(iCell),TimeStepMax)
        dtCell(iCell) = max(dtCell(iCell),TimeStepMin)
!3. dTmin + dtDiff
!        n = int(dtCell(iCell)/TimeStepDiff)
!!!        dtCell(iCell) = TimeStepDiff * float(n)
!        if(n.le.0) then
!!            write(*,*) '@Warning: dtCell<TimeStepMin ! '
!!            n = 1
!            n = int(dtCell(iCell)/TimeStepMin)
!            n = max(1,n)
!
!!!            dtCell(iCell) = TimeStepMin * float(n)
!        else
!            n = nTimeStepDiff * n ! int(dtCell(iCell)/TimeStepMin)
!        end if
!
!
!        if(nNextTimeCell(iCell).eq.nDtGlobal)then
!            CellGlobalStepDone(iCell) = .true.
!        else
!            n = min(nDtGlobal-nNextTimeCell(iCell),n)
!        end if
!        dtCell(iCell) = TimeStepMin * float(n)


!2. 2^n
        n = int(Log(dtCell(iCell)/TimeStepMin + TimeEpsilon)/Log(2.0d0))
        n = max(0,n)
!!!        dtCell(iCell) = TimeStepDiff * float(n)
!!        if(n.le.0) then
!!            write(*,*) '@Warning: dtCell<TimeStepMin ! '
!!            n = 1
!!!            dtCell(iCell) = TimeStepMin * float(n)
!!!        else
!!!            n = nTimeStepDiff * n ! int(dtCell(iCell)/TimeStepMin)
!!        end if

        n = 2**n
        
!        n = int((dtCell(iCell) + TimeEpsilon)/TimeStepMin)
!        n = max(1,n)
!!        n = min(n,nDtGlobal)

        if(nNextTimeCell(iCell).eq.nDtGlobal)then
            CellGlobalStepDone(iCell) = .true.
        else
            n = min(nDtGlobal-nNextTimeCell(iCell),n)
        end if

!        if(nNextTimeCell(iCell).ge.nDtGlobal)then
!            CellGlobalStepDone(iCell) = .true.
!            nNextTimeCell(iCell) = nDtGlobal
!        else
!            n = max(1,n)
!        end if

!        if(nNextTimeCell(iCell).eq.nDtGlobal)then
!            CellGlobalStepDone(iCell) = .true.
!        else
!            n = min(nDtGlobal-nNextTimeCell(iCell),n)
!        end if

        dtCell(iCell) = TimeStepMin * float(n)

!1.
!        n = int(dtCell(iCell)/TimeStepDiff)
!!        dtCell(iCell) = TimeStepDiff * float(n)
!        if(n.le.0) then
!!            write(*,*) '@Warning: dtCell<TimeStepMin ! '
!            n = int(dtCell(iCell)/TimeStepMin)
!            n = max(1,n)
!!            dtCell(iCell) = TimeStepMin * float(n)
!        else
!            n = nTimeStepDiff * n ! int(dtCell(iCell)/TimeStepMin)
!        end if
!
!        if(nNextTimeCell(iCell).eq.nDtGlobal)then
!            CellGlobalStepDone(iCell) = .true.
!        else
!            n = min(nDtGlobal-nNextTimeCell(iCell),n)
!        end if
!        dtCell(iCell) = TimeStepMin * float(n)
!------------------------------------------------------!

!        if(dtGlobal-TimeCell(iCell).lt.TimeStepMin)then
!            CellGlobalStepDone(iCell) = .true.
!        else
!!            dtCell(iCell) = min(dtCell(iCell),TimeGlobal-TimeCell(iCell))
!!            dtCell(iCell) = TimeStepMin * float(n)
!!            dtCell(iCell) = min(dtCell(iCell),dtGlobal-TimeCell(iCell))
!            if(dtCell(iCell).gt.dtGlobal-TimeCell(iCell))then
!                n = int(dtGlobal-TimeCell(iCell)/TimeStepMin)
!                dtCell(iCell) = float(n) * TimeStepMin
!            end if
!        end if

!0.
!        n = int(dtCell(iCell)/TimeStepMin)
!!        n = min(n,nDtGlobal - nNextTimeCell(iCell))
!        dtCell(iCell) = TimeStepMin * float(n)
!        if(n.le.0) then
!            write(*,*) '@Warning: dtCell<TimeStepMin ! '
!            dtCell(iCell) = min(TimeStepMin,dtGlobal-TimeCell(iCell))
!        end if
!
!        nDtGlobal = int(dtGlobal/TimeStepMin)
!
!!!        if(TimeGlobal-TimeCell(iCell).lt.TimeEpsilon)then ! TimeStepMin)then
!        if(nNextTimeCell(iCell).ge.nDtGlobal)then
!!        if(dtGlobal-TimeCell(iCell).lt.TimeEpsilon)then
!            CellGlobalStepDone(iCell) = .true.
!        else
!!            dtCell(iCell) = min(dtCell(iCell),TimeGlobal-TimeCell(iCell))
!            n = min(n,nDtGlobal - nNextTimeCell(iCell))
!            dtCell(iCell) = TimeStepMin * float(n)
!!            dtCell(iCell) = min(dtCell(iCell),dtGlobal-TimeCell(iCell))
!        end if

		dt2Cell(iCell)= 0.5d0 * dtCell(iCell)

        if(AdaptiveMeshRefinement)then
            call CalcAdaptiveMeshRefinement(iCell)
        end if

		NextTimeCell(iCell) = NextTimeCell(iCell) + dtCell(iCell)
!		nNextTimeCell(iCell) = int(TimeCell(iCell)/TimeStepMin) + n
		nNextTimeCell(iCell) = nNextTimeCell(iCell) + n
		dnTimeCell(iCell) = n

!        do numFace = 1,6
!            iFace = meshCells(numFace,iCell)
!            if(iFace.gt.nInternalFaces)then
!                Balans0(1:5) = Balans0(1:5) + fluxFaces(1:5,1,iFace) * dt2Cell(iCell)
!            end if
!        end do
	end do

!	dtCell = minVal(dtCell)
!	dt2Cell = 0.5d0 * dtCell(1)

!	do iCell = 1,nCells
!
!	    if(CellGlobalStepDone(iCell)) cycle
!
!		SubStepCellVars(2:5,iCell) = SubStepCellVars(1,iCell) * SubStepCellVars(2:5,iCell) - &
!                                     sumCellFluxes(2:5,iCell) * dt2Cell(iCell) / meshV(iCell)
!        SubStepCellVars(1,iCell) = SubStepCellVars(1,iCell) - sumCellFluxes(1,iCell) * dt2Cell(iCell) / meshV(iCell)
!        SubStepCellVars(2:5,iCell) = SubStepCellVars(2:5,iCell) / SubStepCellVars(1,iCell)
!
!        Pres = PressureAllVars(SubStepCellVars(1:5,iCell))
!
!!        CharsCellPresent(iCell) = .false.
!    end do

end subroutine CalcCorrectorStep

subroutine CalcAdaptiveMeshRefinement(iCell)
use Data_Mesh
use DataTime, only: dtCell, nNextTimeCell, TimeStepMin
implicit none

    integer,intent(in) :: iCell

    integer :: iN, nCellNeighboursLocal, iCellNeighbour, jCell, n
    real(8) :: dt

    dt = dtCell(iCell) / DbCells(iCell)
    n = nNextTimeCell(iCell)

    nCellNeighboursLocal = nCellNeighbours(iCell)
    do iCellNeighbour = 1,nCellNeighboursLocal
        jCell = CellNeighbours(iCellNeighbour,iCell)
        
        iN = n + int(dt*DbCells(jCell)/TimeStepMin)
        if(iN.ge.nNextTimeCell(jCell)) cycle

        call ReduceTimeStep()
    end do

contains

    subroutine ReduceTimeStep()
    use Data_cfoam, only: Unew,RhoNew,Enew
    use DataCells, only: SubStepCellVars
    use DataFaces, only: CharsNew,CharsOld
    use DataTime, only: NextTimeCell
    implicit none
    
        integer :: i, jNdt, jNtime, jCellFace, jFace, jCellOwner, jSide
        real(8) :: ratio, CellVarsOld(5)
    
        jNdt = int(dtCell(jCell)/TimeStepMin)
        jNtime = nNextTimeCell(jCell)
        ratio = 1.0d0 - float(jNtime - iN) / jNdt ! ratio = 1.0d0 - (nNextTimeCell(jCell) - nNextTimeCell(iCell)) / jNdt
        
        CellVarsOld(1) = RhoNew(iCell)
        CellVarsOld(2:4) = Unew(1:3,iCell)
        CellVarsOld(5) = Enew(iCell)

        ! Cell variables reduction
        do i = 1,5
            SubStepCellVars(i,iCell) = ReduceVar(CellVarsOld(i),SubStepCellVars(i,iCell),ratio)
        end do

        ! Characteristics reduction
        do jCellFace = 1,6
            jFace = meshCells(jCellFace,jCell)
            jCellOwner = meshFaceOwner(jFace)
            if(jCellOwner.eq.jCell)then
                jSide = 1
            else
                jSide = 2
            end if
            do i = 1,5
                CharsNew(i,jSide,jFace) = ReduceVar(CharsOld(i,jSide,jFace),CharsNew(i,jSide,jFace),ratio)
            end do
        end do

        ! TimeCell reduction
!        n = int(dtCell(jCell)/TimeStepMin)
!		NextTimeCell(jCell) = NextTimeCell(jCell) - dtCell(jCell)
!		nNextTimeCell(jCell) = nNextTimeCell(jCell) - n        

        dtCell(jCell) = ratio * dtCell(jCell)
!        n = int(dtCell(jCell)/TimeStepMin)
!
!		NextTimeCell(jCell) = NextTimeCell(jCell) + dtCell(jCell)
!		nNextTimeCell(jCell) = nNextTimeCell(jCell) + n        

		NextTimeCell(jCell) = NextTimeCell(iCell)
		nNextTimeCell(jCell) = nNextTimeCell(iCell)

    end subroutine ReduceTimeStep

    real(8) function ReduceVar(VarOld,VarNew,ratio)
    implicit none
        real(8),intent(in) :: VarOld,VarNew,ratio

        ReduceVar = VarOld + ratio * (VarNew - VarOld)

    end function ReduceVar

end subroutine CalcAdaptiveMeshRefinement
